small project

FA=function(Y,K=2,sigma,epsilon = 1e-3,max_iter=200){
    set.seed(601)
    n=dim(Y)[1]
    p=dim(Y)[2]

    lambda=matrix(c(rep(1,4),numeric(3),numeric(3),rep(1,4)),nrow=7,ncol=2)+matrix(rnorm(p*K,0,sigma),nrow=p,ncol=K)
    phi=diag(rnorm(p,0,1))

    loss=c()
    
    original=matrix(c(rep(1,4),numeric(3),numeric(3),rep(1,4)),nrow=7,ncol=2)%*%t(matrix(c(rep(1,4),numeric(3),numeric(3),rep(1,4)),nrow=7,ncol=2))
    new=lambda%*%t(lambda)
    diff=norm(new-original,type="F")
    
    t=0

    while(diff>epsilon){
        t=t+1
    
        Varx=diag(1,K)-t(lambda)%*%solve(lambda%*%t(lambda)+phi)%*%lambda
        Ex=Y%*%solve(lambda%*%t(lambda)+phi)%*%lambda
        Extx=n*Varx+t(Ex)%*%Ex
    
        lambda_new=t(Y)%*%Ex%*%solve(Extx)
        

        lambda=lambda_new
    
        phi_new=diag(diag(1/n*(t(Y)%*%Y-2*t(Y)%*%Ex%*%t(lambda)+lambda%*%Extx%*%t(lambda))))
        phi=phi_new

        Sigma_y=lambda%*%t(lambda)+phi+epsilon * diag(p) 

        loss<-c(loss,(-n*p/2)*log(2*pi)-(n/2)*log(abs(det(Sigma_y)))-0.5*sum(diag(solve(Sigma_y)%*%t(Y)%*%Y)))
        new=lambda%*%t(lambda)
        diff=norm(new-original,type="F")
        if(t>max_iter){
          break
        }
    }
    
    #print(paste("Times of iteration=",i))
    return(list(lambda=lambda,phi=phi,loss=loss))
    }

n=100
generator_pca=function(n=100){
    Lambda=matrix(c(rep(1,4),numeric(3),numeric(3),rep(1,4)),nrow=7,ncol=2)
    X=mvrnorm(n,numeric(2),diag(1,2))
    W=mvrnorm(n,numeric(7),diag(0.4,7))
    Y=t(Lambda%*%t(X)+t(W))
    return(Y)
    }

set.seed(601)
Y=generator_pca(100)
A1=FA(Y,2,1,1e-3,70)
#set.seed(601)
A2=FA(Y,2,3,1e-3,70)
#set.seed(601)
A3=FA(Y,2,5,1e-3,70)
#set.seed(601)
A4=FA(Y,2,7,1e-3,70)
A5=FA(Y,2,9,1e-3,70)
A6=FA(Y,2,11,1e-3,70)
df1=NULL
df1=data.frame(iteration=c(1:length(A1$loss)),loss=A1$loss)
df2=data.frame(iteration=c(1:length(A2$loss)),loss=A2$loss)
df3=data.frame(iteration=c(1:length(A3$loss)),loss=A3$loss)
df4=data.frame(iteration=c(1:length(A4$loss)),loss=A4$loss)
df5=data.frame(iteration=c(1:length(A5$loss)),loss=A5$loss)
df6=data.frame(iteration=c(1:length(A6$loss)),loss=A6$loss)
plot_ly(data=df1,x=~iteration,y=~loss,type="scatter",mode="lines+markers",name = "Sigma=1")%>%
  add_trace(x = df1[,1], y = df2[,2], 
            type = "scatter", mode = "lines+markers",name = "Sigma=3")%>%
  add_trace(x = df1[,1], y = df3[,2], 
            type = "scatter", mode = "lines+markers",name = "Sigma=5")%>%
  add_trace(x = df1[,1], y = df4[,2], 
            type = "scatter", mode = "lines+markers",name = "Sigma=7")%>%
  add_trace(x = df1[,1], y = df5[,2], 
            type = "scatter", mode = "lines+markers",name = "Sigma=9")%>%
  add_trace(x = df1[,1], y = df6[,2], 
            type = "scatter", mode = "lines+markers",name = "Sigma=11")